
#include <GL/openglut.h>

#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>

typedef struct{
        int Z;
        float r[3];
} ATOMS;
typedef struct{
        int X,Y;
        int type;
} BONDS;

float color_table[128][3];
char periodic_table[][4] = {
     "H",                                                                                 "He",
     "Li","Be",                                                  "B", "C", "N", "O", "F", "Ne",
     "Na","Mg",                                                  "Al","Si","P", "S", "Cl","Ar",
     "K", "Ca","Sc","Ti","V", "Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr",
     "Rb","Sr","Y", "Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I", "Xe",
     "Cs","Ba","La",
         "Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu",
                    "Hf","Ta","W", "Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn",
     "Fr","Ra","Ac",
         "Tb","Pa","U","Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr",
                    "Rf","Ha"};
int getZ(char *symbol){
    int i;
    for(i=0;i<100;i++)
      if (!stricmp(symbol,periodic_table[i])) return i+1;
    return 0;
}

void load_color_table(void){
     FILE *f = fopen("barvy.txt","r");
     char line[512];
     int i,j;
     float r,g,b;
     
     char tmp[8];
     
     for(i=0;i<128;i++) 
      for(j=0;j<3;j++) 
       color_table[i][j]=0.5;
       
     while(!feof(f)){
       fgets(line,512,f);
       sscanf(line,"%s %f %f %f",tmp,&r,&g,&b);
       i = getZ(tmp);
       color_table[i][0]=r;
       color_table[i][1]=g;
       color_table[i][2]=b;        
       printf("Barva %d %f %f %f\n",i,r,g,b);             
     }     
     fclose(f);
}

void load_molfile(char *fn,int *nA,int *nB,ATOMS **A,BONDS **B,float *scale){
     FILE *f = fopen(fn,"r");
     char line[512],tmp[16];
     ATOMS *Atmp;
     BONDS *Btmp;
     int i,j;
     float rmin[3]={1e20,1e20,1e20},rmax[3]={-1e20,-1e20,-1e20};

     for(i=0;i<4;i++)
       fgets(line,512,f);
 
     sscanf(line,"%3d%3d",nA,nB);
 
     *A = (ATOMS *)malloc(*nA*sizeof(ATOMS));       
     *B = (BONDS *)malloc(*nB*sizeof(BONDS));       
     
     for(i=0;i<*nA;i++){
       fgets(line,512,f);
       Atmp = &(*A)[i];
       sscanf(line,"%10f%10f%10f%s",&Atmp->r[0],&Atmp->r[1],&Atmp->r[2],tmp);
       Atmp->Z=getZ(tmp);
     }

     for(i=0;i<*nB;i++){
       fgets(line,512,f);
       Btmp = &(*B)[i];
       sscanf(line,"%3d%3d%3d",&Btmp->X,&Btmp->Y,&Btmp->type);
       Btmp->X--;
       Btmp->Y--;
     }

     for(i=0;i<*nA;i++)for(j=0;j<3;j++){
            if ((*A)[i].r[j]<rmin[j]) rmin[j]=(*A)[i].r[j];        
            if ((*A)[i].r[j]>rmax[j]) rmax[j]=(*A)[i].r[j];        
     }
     
     for(j=0,*scale=0;j<3;j++)
       if (*scale<rmax[j]-rmin[j]) *scale = rmax[j]-rmin[j];
     
     for(i=0;i<*nA;i++)for(j=0;j<3;j++)
            (*A)[i].r[j]=((*A)[i].r[j]-0.5*(rmax[j]+rmin[j]));        
     fclose(f);
}
#define KONSTANTA 0.08
#define RADIUS 1.6

void cross_product(float *c, float *a, float *b){
     c[0] =  a[1]*b[2] - a[2]*b[1];
     c[1] =  a[2]*b[0] - a[0]*b[2];
     c[2] =  a[0]*b[1] - a[1]*b[0];
}

void cross_product_atoms(float *c,ATOMS *A,int a0,int a1, int b0,int b1){
     int k;
     float a[3],b[3];
     for(k=0;k<3;k++){
          a[k]=A[a1].r[k]-A[a0].r[k];            
          b[k]=A[b1].r[k]-A[b0].r[k];            
     }
     cross_product(c,a,b);
}
float dot_product(float *a,float *b){
      return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
}

void normalize_vector(float *x){
     float d;
     int i;
     d = sqrt(dot_product(x,x));
     for(i=0;i<3;i++) x[i]/=d;
}

void display_model(int nA,int nB,ATOMS *A,BONDS *B, float scale){
     int i,j,k,bf,b1[4],b2[4];
     float r[3],d,ut[4][3],u[3],l[3];
     for(i=0;i<nA;i++){
           glPushMatrix();
           glTranslatef(A[i].r[0],A[i].r[1],A[i].r[2]);
           printf("Atom %d %f %f %f\n",A[i].Z,A[i].r[0],A[i].r[1],A[i].r[2]);
           glColor3f(color_table[A[i].Z][0],color_table[A[i].Z][1],color_table[A[i].Z][2]);
           printf("Color %f %f %f\n",color_table[A[i].Z][0],color_table[A[i].Z][1],color_table[A[i].Z][2]);
           glutSolidSphere(KONSTANTA*3,20,20);
           glPopMatrix();            
     }
     glColor3f(0.5,.5,.5);
     for(i=0;i<nB;i++){
           glPushMatrix();
           glTranslatef(A[B[i].X].r[0],A[B[i].X].r[1],A[B[i].X].r[2]);
           for(j=0;j<3;j++) r[j]=A[B[i].Y].r[j]-A[B[i].X].r[j];
           d = sqrt(dot_product(r,r));
           switch(B[i].type){
                default:             
                case 1:
                   glRotatef(180,r[0],r[1],r[2]+d);
                   glutSolidCylinder(KONSTANTA,d,20,20);
                   break;
                case 3:
                   glRotatef(180,r[0],r[1],r[2]+d);
                   glTranslatef(RADIUS*KONSTANTA,0,0);
                   glutSolidCylinder(KONSTANTA,d,20,20);
                   glTranslatef(-RADIUS*(1.5)*KONSTANTA,RADIUS*0.866025404*KONSTANTA,0);
                   glutSolidCylinder(KONSTANTA,d,20,20);
                   glTranslatef(0,-2*RADIUS*0.866025404*KONSTANTA,0);
                   glutSolidCylinder(KONSTANTA,d,20,20);                   
                   break;                    
               case 2:
                    printf("Bond %d\n",i);
                    for(j=bf=0;j<nB;j++)if(i!=j){
                         if (B[i].X==B[j].X){
                            b1[bf]=B[i].X;
                            b2[bf++]=B[j].Y;
                            printf("bf%d(%d-%d)\n",bf,b1[bf-1],b2[bf-1]);
                         }
                         if (B[i].X==B[j].Y){
                            b1[bf]=B[i].X;
                            b2[bf++]=B[j].X;
                            printf("bf%d(%d-%d)\n",bf,b1[bf-1],b2[bf-1]);
                         }
                         if (B[i].Y==B[j].X){
                            b1[bf]=B[i].Y;
                            b2[bf++]=B[j].Y;
                            printf("bf%d(%d-%d)\n",bf,b1[bf-1],b2[bf-1]);
                         }
                         if (B[i].Y==B[j].Y){
                            b1[bf]=B[i].Y;
                            b2[bf++]=B[j].X;
                            printf("bf%d(%d-%d)\n",bf,b1[bf-1],b2[bf-1]);
                         }
                    }
                    for(j=0;j<bf;j++){
                       cross_product_atoms(ut[j],A,B[i].X,B[i].Y,b1[j],b2[j]);
                       normalize_vector(ut[j]);               
                    }
                    if (bf) for(k=0;k<3;k++) u[k]=ut[0][k];
                    for(j=1;j<bf;j++){
                         if (dot_product(u,ut[j])>0){
                              for(k=0;k<3;k++) u[k]+=ut[j][k];
                         }else{
                              for(k=0;k<3;k++) u[k]-=ut[j][k];
                         }                              
                    }
                    if(bf){
                         cross_product(l,u,r);
                         normalize_vector(l);
                         glPushMatrix();
                         glTranslatef(RADIUS*KONSTANTA*l[0],RADIUS*KONSTANTA*l[1],RADIUS*KONSTANTA*l[2]);
                         glRotatef(180,r[0],r[1],r[2]+d);
                         glutSolidCylinder(KONSTANTA,d,20,20);
                         glPopMatrix();  
                         glPushMatrix();
                         glTranslatef(-RADIUS*KONSTANTA*l[0],-RADIUS*KONSTANTA*l[1],-RADIUS*KONSTANTA*l[2]);
                         glRotatef(180,r[0],r[1],r[2]+d);
                         glutSolidCylinder(KONSTANTA,d,20,20);
                         glPopMatrix();  
                    } else {
                         glPushMatrix();
                         glTranslatef(RADIUS*KONSTANTA,0,0);
                         glRotatef(180,r[0],r[1],r[2]+d);
                         glutSolidCylinder(KONSTANTA,d,20,20);
                         glPopMatrix();  
                         glPushMatrix();
                         glTranslatef(-RADIUS*KONSTANTA,0,0);
                         glRotatef(180,r[0],r[1],r[2]+d);
                         glutSolidCylinder(KONSTANTA,d,20,20);
                         glPopMatrix();    
                    }
                   break;     
           }
           glPopMatrix();            
     }

}

void take_a_picture(void){
     unsigned char *data;
     FILE *f=fopen("obrazek.bmp","wb");
     char header[54]={     
          66, 77, 54, -32, 46, 0, 0, 0, 
          0, 0, 54, 0, 0, 0, 40, 0,
          0, 0, 0, 5, 0, 0, 32, 3,
          0, 0, 1, 0, 24, 0, 0, 0, 0, 0, 0, -32, 46, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0
     };     
     int w,h,*i;    
     w = glutGet(GLUT_WINDOW_WIDTH);
     h = glutGet(GLUT_WINDOW_HEIGHT);
/*
     i = (void *)(&header[17]);
     *i = w;
     i = (void *)(&header[21]);
     *i = h;
*/
     data = (unsigned char *)malloc(w*h*3);
     glReadBuffer(GL_BACK); 
     glPixelZoom(1.0f,1.0f);
     glReadPixels(0,0,w,h,GL_RGB,GL_UNSIGNED_BYTE, data); 
     fwrite(header,1,sizeof(header),f);
     fwrite(data,1,w*h*3,f);
     fclose(f);
}


/* data */
int nA,nB;
ATOMS *A;
BONDS *B;
float scale;

void resize(int width, int height)
{
    const float ar = (float) width / (float) height;
    glViewport(0, 0, width, height);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    if( ar > 1 )
        glFrustum( -ar, ar, -1.0, 1.0, 2.0, 100.0 );
    else
        glFrustum( -1.0, 1.0, -1/ar, 1/ar, 2.0, 100.0 );
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity() ;
}

int list;
void display(void)
{
    const double t = glutGet(GLUT_ELAPSED_TIME) / 1000.0;

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glLoadIdentity();
    glTranslatef(0,0,-6);
    if (!glIsList(list)){
         glNewList(list=glGenLists(1),GL_COMPILE);
         glScalef(4/scale,4/scale,4/scale);
         display_model(nA,nB,A,B,scale);         
         glEndList();                
    }
//    glRotatef(10*t,1,1,1);
    glCallList(list);
    glutSwapBuffers();
}


void key(unsigned char key, int x, int y)
{
    switch (key){
    case 27 :
    case 'Q':
    case 'q':
        exit(0);
        break;
    case 's':
    case 'S':
         take_a_picture();
         break;    
    default:
        break;
    }
    glutPostRedisplay();
}


void idle(void){
    glutPostRedisplay();
}

const GLfloat light_ambient[]  = { 0.0f, 0.0f, 0.0f, 1.0f };
const GLfloat light_diffuse[]  = { 1.0f, 1.0f, 1.0f, 1.0f };
const GLfloat light_specular[] = { 1.0f, 1.0f, 1.0f, 1.0f };
const GLfloat light_position[] = { 2.0f, 5.0f, 5.0f, 0.0f };

const GLfloat mat_ambient[]    = { 0.7f, 0.7f, 0.7f, 1.0f };
const GLfloat mat_diffuse[]    = { 0.8f, 0.8f, 0.8f, 1.0f };
const GLfloat mat_specular[]   = { 1.0f, 1.0f, 1.0f, 1.0f };
const GLfloat high_shininess[] = { 100.0f };

int main(int argc, char **argv)
{
    glutInitWindowSize(640,480);
    glutInitWindowPosition(40,40);
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);

    glutCreateWindow(argv[1]);

    glutReshapeFunc(resize);
    glutDisplayFunc(display);
    glutKeyboardFunc(key);
    glutIdleFunc(idle);

    glClearColor(1,1,1,1);

    glEnable(GL_DEPTH_TEST);
    glDepthFunc(GL_LESS);

    glEnable(GL_LIGHTING);
    glEnable(GL_LIGHT0);
    glEnable(GL_NORMALIZE);
    glEnable(GL_COLOR_MATERIAL);

    glLightfv(GL_LIGHT0, GL_AMBIENT,  light_ambient);
    glLightfv(GL_LIGHT0, GL_DIFFUSE,  light_diffuse);
    glLightfv(GL_LIGHT0, GL_SPECULAR, light_specular);
    glLightfv(GL_LIGHT0, GL_POSITION, light_position);

    glMaterialfv(GL_FRONT, GL_AMBIENT,   mat_ambient);
    glMaterialfv(GL_FRONT, GL_DIFFUSE,   mat_diffuse);
    glMaterialfv(GL_FRONT, GL_SPECULAR,  mat_specular);
    glMaterialfv(GL_FRONT, GL_SHININESS, high_shininess);

    glHint(GL_PERSPECTIVE_CORRECTION_HINT,GL_NICEST);
    glEnable (GL_POLYGON_SMOOTH);
    glHint(GL_POLYGON_SMOOTH_HINT,GL_NICEST);
    

glEnable(GL_BLEND);
glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

glHint(GL_POINT_SMOOTH, GL_NICEST);
glHint(GL_LINE_SMOOTH, GL_NICEST);
glHint(GL_POLYGON_SMOOTH, GL_NICEST);

glEnable(GL_POINT_SMOOTH);
glEnable(GL_LINE_SMOOTH);
glEnable(GL_POLYGON_SMOOTH);


    load_color_table();
    load_molfile(argv[1],&nA,&nB,&A,&B,&scale);
    glutFullScreen();
    glutMainLoop();
    return EXIT_SUCCESS;
}
